#!/usr/bin/python
"""
Processes associations, assumes that their features have been processed, associations tsv file must be in the format of
feature1\tfeature2\tpvalue\timportance\tcorrelation 
"""
import sys
import db_util
import os
import time
import parse_features_rfex

myhost = db_util.getDBHost() 
mydb = db_util.getDBSchema() 
myuser = db_util.getDBUser() 
mypw = db_util.getDBPassword() 

args = sys.argv
pvalue_cutoff = 1.0
if (len(args) < 4):
	print "Usage is py2.6 parse_associations_rfex.py data_matrix.tsv data_associations.tsv dataset_label [pvalue_cutoff]"
	sys.exit(1)

matrixfile = args[1]
associationsfile = args[2]
dataset_label = args[3]
associations_table = mydb + "." + dataset_label + "_networks"
features_table = mydb + "." + dataset_label + "_features"
if (len(args) >= 5):
	pvalue_cutoff = float(args[4])

print "Begin processing associations %s Applying pvalue cutoff %f" %(time.ctime(), pvalue_cutoff)
associations_in = open(associationsfile,'r')
matrix_in = open(matrixfile,'r')
featureId = 0
features_hash = {}
parse_features_rfex.process_feature_matrix(matrixfile)
#print "number of features in parse:" + str(len(parse_features_rfex.features_hash))
path = sys.argv[1].rsplit("/", 1)[0]
if (os.path.exists(path + "/GEXP_interestingness.tsv")):
	parse_features_rfex.process_gexp_interest_score(path + "/GEXP_interestingness.tsv")
#print "number of feature gene interesting:" + str(len(parse_features_rfex.gene_interesting_hash))

fshout = open('./results/load_sql_associations_' + dataset_label + '.sh','w')
tsvout = open('./results/' + dataset_label + '_links.tsv','w')

def isUnmappedAssociation(f1alias, f2alias):
	"""
	Skipping N:METH:cg16793061:chr8:72436538:: B:SAMP:I(C1|mRNA_CC_3)::::
	"""	
	f1data = f1alias.split(":")
	f1source = f1data[1]
	f1chr = f1data[3]
	f2data = f2alias.split(":")
	f2source = f2data[1]
	f2chr = f2data[3]
	
	if (f1chr == "" and (f1source != "CLIN" and f1source != "SAMP")):
		if (f2source != "CLIN" and f2source != "SAMP"):
			return True
	if (f2chr == "" and (f2source != "CLIN" and f2source != "SAMP")):
		if (f1source != "CLIN" and f1source != "SAMP"):
			return True
	return False	

lc = 0
unMapped = 0
pvalueCutCount = 0
lines = associations_in.readlines()
associations_in.close()
for line in lines:
	lc = lc + 1
	columns = line.strip().split('\t')
	f1alias = columns[0]
	f2alias = columns[1]
	pvalue = columns[2]
	importance = columns[3]
	correlation = columns[4]
	if (isUnmappedAssociation(f1alias, f2alias)):
		#print "Skipping %s %s" %(f1alias, f2alias)
		unMapped += 1
		continue	
	if (not parse_features_rfex.getFeatureId(f1alias) or not parse_features_rfex.getFeatureId(f2alias)):
		print "feature1 %s or feature2 %s not in original feature matrix" %(f1alias, f2alias)
		continue
	f1genescore = ""
	if (parse_features_rfex.getGeneInterestScore(f1alias) != None):
		f1genescore = parse_features_rfex.getGeneInterestScore(f1alias)	
	f2genescore = ""
	if (parse_features_rfex.getGeneInterestScore(f2alias) != None):
		f2genescore = parse_features_rfex.getGeneInterestScore(f2alias)
	f2data = f2alias.split(':')
	if len(f2data) > 4:
		f2data[3] = f2data[3][3:]
	if (float(pvalue) <= pvalue_cutoff): 
		tsvout.write(str(lc) + "\t" + f1alias + "\t" + f2alias + "\t" + pvalue + "\t" + importance + "\t" + correlation + "\t" + str(parse_features_rfex.getFeatureId(f1alias)) + "\t" + str(parse_features_rfex.getFeatureId(f2alias)) + "\t" + "\t".join(f2data) + "\t" + f1genescore + "\t" + f2genescore + "\n")
	else:
		pvalueCutCount += 1

fshout.write("#!/bin/bash\n")
fshout.write("mysql --user=%s --password=%s --database=%s<<EOFMYSQL\n" %(myuser, mypw, mydb))
fshout.write("load data local infile './results/" + dataset_label + "_links.tsv" + "' replace INTO TABLE " + associations_table + " fields terminated by '\\t' LINES TERMINATED BY '\\n';")
fshout.write("\nEOFMYSQL\n")

tsvout.close()
fshout.close()
print "%i associations filtered because of > than pvalue %f, unMapped %i" %(pvalueCutCount, pvalue_cutoff, unMapped)
print "Begin bulk upload %s os.system sh %s" %(time.ctime(), fshout.name)
os.system("sh " + fshout.name)
print "Done processing associations %s" %(time.ctime())

